#! /cygdrive/c/PROGRA~1/R/R-2.12.1/bin/R

library(RODBC)
library(DESeq)
library(biomaRt)
library(sqldf)
conn<- odbcConnect(dsn= 'pgVitelleschi')

sqlQuery(conn, "SELECT cross_tab('SELECT gene_id, source, \"count\" FROM htseq_count', 'tmp_htseq')")
countsTable<- sqlQuery(conn, 'select * from tmp_htseq')
sqlQuery(conn, 'drop table tmp_htseq')
genes<- countsTable[,1]
countsTable<- countsTable[,2:ncol(countsTable)]
rownames(countsTable)<- genes
countsTable[is.na(countsTable)]<- 0
head(countsTable)

# -----------------------------------------------------------------------------
#                                  edgeR: CTRL vs LPS
# -----------------------------------------------------------------------------

library(edgeR)
## See edgeR vignette edger.pdf page 51 "Case study: Oral carcinomas vs matched normal tissue"
targets<- as.data.frame(cbind(
    treatment= c('ctrl', 'lps', 'ctrl', 'lps'),
    mphage= c('am', 'am', 'bmdm', 'bmdm')
))
row.names(targets)<- paste(targets$mphage, targets$treatment, sep= '_')
targets
#           treatment mphage
# am_ctrl        ctrl     am
# am_lps          lps     am
# bmdm_ctrl      ctrl   bmdm
# bmdm_lps        lps   bmdm

d<- countsTable
colnames(d) <- row.names(targets)
d<- DGEList(counts = d, group = targets$treatment)
d<- calcNormFactors(d)
d<- d[rowSums(d$counts) > 9, ] ## Remove genes with tot number of reads < 10

plotMDS.dge(d, main="MDS Plot for CTRL-LPS data")
## Note: LPS tratement doesn't seem to have more effect than tissue of origin (AM vs BMDM). Is this suspicious?
savePlot('M:/Documents/LabBook/LabBook_Figures/20110712_MDSplot_edger.emf', 'emf')
cor(countsTable)
#                    20110202_am_ctrl 20110202_am_lps 20110202_bmdm_ctrl 20110202_bmdm_lps
# 20110202_am_ctrl             1.0000          0.8540             0.5882            0.6419
# 20110202_am_lps              0.8540          1.0000             0.4829            0.6162
# 20110202_bmdm_ctrl           0.5882          0.4829             1.0000            0.8989
# 20110202_bmdm_lps            0.6419          0.6162             0.8989            1.0000

design <- model.matrix(~ mphage + treatment, data = targets)
#           (Intercept) mphagebmdm treatmentlps
# am_ctrl             1          0            0
# am_lps              1          0            1
# bmdm_ctrl           1          1            0
# bmdm_lps            1          1            1

d<- estimateCRDisp(d, design)
d$CR.common.dispersion       ## 0.14: commnon dispersion
sqrt(d$CR.common.dispersion) ## 0.37: coefficient of biological variation

glmfit <- glmFit(d, design, dispersion = d$CR.common.dispersion)
lrt.d <- glmLRT(d, glmfit)
options(digits = 4)
degenes<- topTags(lrt.d, n= nrow(lrt.d$table))$table
head(degenes)

countsTable[rownames(countsTable) == 'ENSSSCG00000001404', ] ## This is TNFA

## Export to file and to postgres
write.table(degenes, 'F:/data/20110711_rnaseq-am-bmdm/lps_vs_ctrl.edger.txt', row.names=TRUE, col.names=TRUE, sep= '\t', quote= FALSE)
sqlQuery(conn, "SELECT read_table($$ file: 'F:/data/20110711_rnaseq-am-bmdm/lps_vs_ctrl.edger.txt', table: 'rnaseq_de.edger_toptags', header: ['id', 'logconc', 'logfc', 'pvalue', 'fdr'], skip: 1 $$)")
sqlQuery(conn, "comment on table rnaseq_de.edger_toptags is 'Output from edgeR comparing lps vs ctrl in rnaseq libs AM and BMDM using (paired samples). logfc > 0 means overexpressed in LPS. See rnaseq-am-bmdm/deseq_edger.R' ")

summary(decideTestsDGE(lrt.d))
graphics.off()
windows(height= 12/2.54, width = 12/2.54)
par(cex= 0.85)
plotSmear(lrt.d, de.tags= rownames(degenes[degenes$FDR < 0.05,]), main= 'Genes expression in macrophages: LPS vs CTRL')
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000001404'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000001404'], labels= 'TNFA', cex= 0.85)
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000008088'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000008088'], labels= 'IL1B', cex= 0.85)
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000007585'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000007585'], labels= 'ACTB', cex= 0.85)
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000008090'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000008090'], labels= 'IL1A', cex= 0.85)
savePlot('M:/Documents/LabBook/LabBook_Figures/20110712_expression_edger.tiff', 'tiff')

# -----------------------------------------------------------------------------
#                           goseq
# -----------------------------------------------------------------------------
library(goseq)
genes_sorted<- degenes[order(rownames(degenes)),]
head(genes_sorted)
genes<- as.integer(genes_sorted$FDR < 0.05)
names(genes)<- row.names(genes_sorted)
table(genes)
# genes
#     0     1 
# 10446   767
genes[1:10]

ensembl<-  useMart("ensembl")
mart<- useDataset("sscrofa_gene_ensembl", ensembl)
cdna<- getBM(
    attributes= c("cdna", "ensembl_transcript_id", "ensembl_gene_id"), 
    filters= "ensembl_gene_id",
    value= rownames(degenes),
    mart= mart
)
cdna$cdna_length<- nchar(cdna$cdna)
gene_length<- aggregate(cdna$cdna_length, by= list(cdna$ensembl_gene_id), median)
gene_ids<- gene_length$Group.1
gene_length<- gene_length$x
gene_length[gene_length > 25000] <- NA ## There is some problem with pnull if this is not done (only 1 NA generated)
names(gene_length)<- gene_ids

pwf <- nullp(genes, bias.data= gene_length)

go_genes<- getBM(
    attributes= c("ensembl_gene_id", "name_1006"),
    filters= "ensembl_gene_id",
    value= rownames(degenes),
    mart= mart
)
names(go_genes)<- c("ensembl_gene_id", "go_term_name")
go_genes<- go_genes[go_genes$go_term_name != '', ]
head(go_genes)

GO.wall <- goseq(pwf, gene2cat= go_genes)

## --- Check length bias (none)
GO.hyper<- goseq(pwf, gene2cat= go_genes, method= 'Hypergeometric')

plot(log10(GO.wall[, 2]), log10(GO.hyper[match(GO.hyper[, 1],
    GO.wall[, 1]), 2]), xlab = "log10(Wallenius p-values)", ylab = "log10(Hypergeometric)", xlim = c(-3, 0), ylim = c(-3, 0))
abline(0, 1, col = 3, lty = 2)
## ---

GO.wall$over_represented_fdr<- p.adjust(GO.wall$over_represented_pvalue, method= 'fdr')
GO.wall$under_represented_fdr<- p.adjust(GO.wall$under_represented_pvalue, method= 'fdr')
sqldf('select * from "GO.wall" where under_represented_fdr < 0.05 or over_represented_fdr < 0.05')

## Annotate genes
annotation<- getBM(
    attributes= c("ensembl_gene_id", "external_gene_id", "wikigene_name", "description", "wikigene_description"),
    filters= "ensembl_gene_id",
    value= rownames(degenes),
    mart= mart
)
write.table(annotation, file= 'F:/data/20110711_rnaseq-am-bmdm/gene_annotation.txt', row.names=FALSE, col.names=TRUE, sep= '\t', quote= FALSE)
sqlQuery(conn, "select read_table($$ file: 'F:/data/20110711_rnaseq-am-bmdm/gene_annotation.txt', table: 'rnaseq_de.gene_annotation', header: True, overwrite: True $$) ")
sqlQuery(conn, "COMMENT ON TABLE rnaseq_de.gene_annotation IS 'Annotation of genes in edger_toptags extracted from biomart. See rnaseq-am-bmdm/deseq_edger.R' ")

go_annotation<- sqldf('select distinct 
      annotation.ensembl_gene_id, wikigene_name, description, category, over_represented_fdr, under_represented_fdr
      from "GO.wall" inner join go_genes on category = go_term_name
                     inner join annotation on annotation.ensembl_gene_id = go_genes.ensembl_gene_id
      -- where under_represented_fdr < 0.05 or over_represented_fdr < 0.05
      order by wikigene_name')
write.table(go_annotation, file= 'F:/data/20110711_rnaseq-am-bmdm/go_annotation.txt', row.names=FALSE, col.names=TRUE, sep= '\t', quote= FALSE)
sqlQuery(conn, "select read_table($$ file: 'F:/data/20110711_rnaseq-am-bmdm/go_annotation.txt', table: 'rnaseq_de.goseq_annotation', header: True $$) ")
sqlQuery(conn, "COMMENT ON TABLE rnaseq_de.goseq_annotation IS 'Genes in RNAseq AM and BMDM with GO annotation and over/under representation produced by goseq. See rnaseq-am-bmdm/deseq_edger.R' ")

# -----------------------------------------------------------------------------
#                          edgeR: Tissue D.E. AM vs BMDM
#
# This analysis not used as probably the difference in library prep leads to
# a large number (3406 out of 11213) of DE genes, which is suspicious.
# -----------------------------------------------------------------------------

library(edgeR)
## See edgeR vignette edger.pdf page 51 "Case study: Oral carcinomas vs matched normal tissue"
targets<- as.data.frame(cbind(
    treatment= c('ctrl', 'lps', 'ctrl', 'lps'),
    mphage= c('am', 'am', 'bmdm', 'bmdm')
))
row.names(targets)<- paste(targets$mphage, targets$treatment, sep= '_')
targets
#           treatment mphage
# am_ctrl        ctrl     am
# am_lps          lps     am
# bmdm_ctrl      ctrl   bmdm
# bmdm_lps        lps   bmdm

d<- countsTable
colnames(d) <- row.names(targets)
d<- DGEList(counts = d, group = targets$treatment)
d<- calcNormFactors(d)
d<- d[rowSums(d$counts) > 9, ] ## Remove genes with tot number of reads < 10

design <- model.matrix(~ mphage + treatment, data = targets)
d<- estimateCRDisp(d, design)
d$CR.common.dispersion       ## 0.14: commnon dispersion
sqrt(d$CR.common.dispersion) ## 0.37: coefficient of biological variation

glmfit <- glmFit(d, design, dispersion = d$CR.common.dispersion)
lrt.d <- glmLRT(d, glmfit, coef= 2)
options(digits = 4)
degenes<- topTags(lrt.d, n= nrow(lrt.d$table))$table
head(degenes)
dim(degenes[degenes$FDR < 0.05, ])
dim(degenes)
hist(degenes$FDR)

countsTable[rownames(countsTable) == 'ENSSSCG00000017296', ]

summary(decideTestsDGE(lrt.d))
windows(height= 12/2.54, width = 12/2.54)
par(cex= 0.85)
## +ve logfc: Overexpressed in BMDM
plotSmear(lrt.d, de.tags= rownames(degenes[degenes$FDR < 0.05,]), main= 'Genes expression in macrophages: AM vs BMDM')
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000001404'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000001404'], labels= 'TNFA', cex= 0.85)
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000008088'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000008088'], labels= 'IL1B', cex= 0.85)
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000007585'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000007585'], labels= 'ACTB', cex= 0.85)
text(x= degenes$logConc[rownames(degenes) == 'ENSSSCG00000008090'], y= degenes$logFC[rownames(degenes) == 'ENSSSCG00000008090'], labels= 'IL1A', cex= 0.85)


# -----------------------------------------------------------------------------
#                           GOseq on AM vs BMDM
# -----------------------------------------------------------------------------

library(goseq)
genes_sorted<- degenes[order(rownames(degenes)),]
head(genes_sorted)
genes<- as.integer(genes_sorted$FDR < 0.05)
names(genes)<- row.names(genes_sorted)
table(genes)

genes[1:10]

ensembl<-  useMart("ensembl")
mart<- useDataset("sscrofa_gene_ensembl", ensembl)
cdna<- getBM(
    attributes= c("ensembl_transcript_id", "ensembl_gene_id", "cdna"), 
    filters= "ensembl_gene_id",
    value= rownames(degenes),
    mart= mart
)
cdna$cdna_length<- nchar(cdna$cdna)
cdna<- cdna[, c("ensembl_transcript_id", "ensembl_gene_id", "cdna_length", "cdna")]
write.table(cdna, file= 'D:/Tritume/rnaseq-am-bmdm/cdna_sscrofa_gene_ensembl.txt', row.names= FALSE, col.names= TRUE, sep= '\t', quote= FALSE)
system('gzip D:/Tritume/rnaseq-am-bmdm/cdna_sscrofa_gene_ensembl.txt')
gene_length<- aggregate(cdna$cdna_length, by= list(cdna$ensembl_gene_id), median)
gene_ids<- gene_length$Group.1
gene_length<- gene_length$x
gene_length[gene_length > 25000] <- NA ## There is some problem with pnull if this is not done (only 1 NA generated)
names(gene_length)<- gene_ids

pwf <- nullp(genes, bias.data= gene_length)

go_genes<- getBM(
    attributes= c("ensembl_gene_id", "name_1006"),
    filters= "ensembl_gene_id",
    value= rownames(degenes),
    mart= mart
)
names(go_genes)<- c("ensembl_gene_id", "go_term_name")
go_genes<- go_genes[go_genes$go_term_name != '', ]
head(go_genes)

GO.wall <- goseq(pwf, gene2cat= go_genes)

## --- Check length bias (none)
GO.hyper<- goseq(pwf, gene2cat= go_genes, method= 'Hypergeometric')

plot(log10(GO.wall[, 2]), log10(GO.hyper[match(GO.hyper[, 1],
    GO.wall[, 1]), 2]), xlab = "log10(Wallenius p-values)", ylab = "log10(Hypergeometric)", xlim = c(-3, 0), ylim = c(-3, 0))
abline(0, 1, col = 3, lty = 2)
## ---

GO.wall$over_represented_fdr<- p.adjust(GO.wall$over_represented_pvalue, method= 'fdr')
GO.wall$under_represented_fdr<- p.adjust(GO.wall$under_represented_pvalue, method= 'fdr')
sqldf('select * from "GO.wall" where under_represented_fdr < 0.05 or over_represented_fdr < 0.05')


# -----------------------------------------------------------------------------
#                             DESeq (Not used)
# -----------------------------------------------------------------------------

conds<- c('ctrl', 'lps', 'ctrl', 'lps')
#conds<- c('am', 'am', 'bmdm', 'bmdm')
cds <- newCountDataSet( countsTable, conds )
cds <- estimateSizeFactors( cds )
cds <- estimateVarianceFunctions( cds, method= 'normal' )

res <- nbinomTest( cds, "ctrl", "lps")
#res <- nbinomTest( cds, "am", "bmdm")
length(res$padj[res$padj < 0.05])
head(res)
plot( res$baseMean, res$log2FoldChange, log="x", pch=20, cex=.1, col = ifelse( res$padj < .1, "red", "black" ) )

ensembl<-  useMart("ensembl")
mart<- useDataset("sscrofa_gene_ensembl", ensembl)
genes<- getBM(
    attributes= c("ensembl_gene_id", "external_gene_id", "wikigene_name", "description"),
    filters= "ensembl_gene_id",
    value= res$id[res$padj < 0.05],
    mart= mart
)
write.table(genes, file= 'F:/data/20110711_rnaseq-am-bmdm/am_bmdm.deseq.txt', row.names=FALSE, col.names=TRUE, sep= '\t', quote= FALSE)
genes[1:100,]
listAttributes(mart)


# ------------------------------[ Tritume ]------------------------------------

plot(countsTable[, c(1,2)])
cor(countsTable)

plot(countsTable[, c(3,4)])
cor(countsTable[, c(3,4)])

plot(countsTable[, c(1,3)])
cor(countsTable[, c(1,3)])

d<- countsTable
head(d)
group<- c('ctrl', 'lps', 'ctrl', 'lps')
d<- DGEList(counts = d, group = group)
d

d <- estimateCommonDisp(d)
d
d$samples$lib.size
d$common.lib.size
colSums(d$pseudo.alt)
d$common.dispersion
sqrt(d$common.dispersion)

de.com <- exactTest(d)
topTags(de.com)
gene_id<- gene_length$ensembl_gene_id
glength<- gene_length$end_position - gene_length$start_position
names(glength)<- gene_id
glength[1:10]

go_genes<- getBM(
    attributes= c("ensembl_gene_id", "name_1006"),
    filters= "ensembl_gene_id",
    value= rownames(degenes)[1:10],
    mart= mart
)
names(go_genes)<- c("ensembl_gene_id", "go_term_name")
go_genes<- go_genes[go_genes$go_term_name != '', ]
head(go_genes)


sg<- supportedGenomes()
sg[sg$db == 'susScr2', ]
pwf <- nullp(genes, "susScr2", "ensGene")

genes<- getBM(
    attributes= c("ensembl_gene_id", "external_gene_id", "hgnc_id", "hgnc_symbol", "wikigene_name", "name_1006", "description"),
    filters= "ensembl_gene_id",
    value= rownames(degenes),
    mart= mart
)
names(genes)
listAttributes(mart)
genes[1:100,]
dim(genes)

writeClipboard(rownames(degenes[degenes$FDR < 0.05 & degenes$logFC > 0,]))

genes<- ifelse(degenes$FDR<0.05, 1, 0)
genes[1:10]


genes<- c('ENSG00000252775', 'ENSG00000207459', 'ENSG00000252899', 'ENSG00000201298', 'ENSG00000222266', 'ENSG00000222924', 'ENSG00000212460', 'ENSG00000199880', 'ENSG00000252533', 'ENSG00000200378')
gene_lengths <- getlength(names(genes), "hg18", "ensGene")